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Abstract 

A  curvilinear  version  of  the  nearshore  circulation  model  SHORECIRC  is  developed  based  on  the  quasi-3D  nearshore 
circulation  equations  derived  by  Putrevu  and  Svendsen  [Eur.  J.  Mech.  18  (1999)  409-427].  We  use  a  generalized  coordinate 
transformation  and  re-derive  the  equations  in  tensor-invariant  forms.  The  contravariant  component  technique  is  used  to  simplify 
both  the  transformed  equations  and  boundary  conditions.  A  high-order  finite-difference  scheme  with  a  staggered  grid  in  the 
image  domain  is  adopted  for  the  numerical  model.  Very  good  convergence  rates  with  both  grid  refinement  and  time  refinement 
are  obtained  in  a  simple  convergence  test.  The  model  is  then  applied  to  four  cases  involving  either  a  non-orthogonal  quadrilateral 
grid  or  a  generalized  curvilinear  grid.  The  versatility  of  the  curvilinear  model  in  dealing  with  curved  shorelines,  nearshore 
breakwaters  and  other  complicated  geometries  is  demonstrated  in  the  test  cases.  The  accuracy  of  the  model  is  shown  in  the  paper 
through  model/data  comparisons  in  two  of  the  case  studies. 

©  2003  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

The  performance  of  a  numerical  model  for  wave- 
induced  nearshore  circulation  is  important  for  predict¬ 
ing  sediment  and  pollutant  transport  in  coastal  regions. 
The  model  performance  mainly  depends  on  the  under¬ 
standing  of  nearshore  phenomena  as  well  as  numerical 
techniques  used  in  the  model.  Over  the  last  decades, 
significant  progress  has  been  made  in  our  understand¬ 
ing  of  wave-generated  phenomena  such  as  wave  set-up, 
wave  breaking,  undertow,  cross-shore  and  longshore 
currents  and  their  stability,  turbulence  and  mixing,  and 
the  generation  of  long-wave  phenomena.  Based  on  this 
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understanding,  a  variety  of  numerical  models  have 
been  developed  and  used  for  modeling  of  nearshore 
phenomena. 

Recently,  the  three-dimensional  dispersion  of  mo¬ 
mentum  in  wave-induced  nearshore  currents  was  dis¬ 
cussed  by  Svendsen  and  Putrevu  (1994).  They  found 
that  the  vertical  nonuniformity  of  the  currents  leads  to 
a  mixing-like  term  in  the  depth-integrated  alongshore 
momentum  equation,  which  is  analogous  to  the  shear- 
dispersion  mechanism  found  by  Taylor  (1953,  1954). 
The  lateral  mixing  caused  by  the  shear-dispersion 
mechanism  is  an  order  of  magnitude  larger  than  the 
turbulent  lateral  mixing  and  is  thus  considered  to  be  a 
major  contributor  to  the  total  lateral  mixing  in  the 
nearshore  region.  Smith  (1997)  also  presented  a  rather 
general  derivation  of  the  shear  dispersion  mechanism 
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for  the  case  with  no  short- wave-induced  volume  flux. 
Putrevu  and  Svendsen  (1999,  PS99  hereafter)  extend¬ 
ed  the  results  of  Svendsen  and  Putrevu  (1994)  to  the 
general  case  of  unsteady  circulations  induced  by  wave 
breaking  over  an  arbitrary  bottom  topography.  For  a 
zero  wave-induced  volume  flux,  the  newly  derived 
equations  are  similar  to  those  of  Smith.  A  quasi-3D 
numerical  model  named  SHORECIRC  (Svendsen  et 
al.,  2000)  has  been  developed  based  on  the  nearshore 
circulation  equations.  It  is  a  2D  horizontal  model 
which  incorporates  the  effect  of  the  vertical  structure 
of  horizontal  flows.  A  semi-analytical  solution  is  used 
for  the  3D  current  profiles  in  combination  with  a 
numerical  solution  for  the  depth- integrated  2D  hori¬ 
zontal  equations.  Several  applications  of  the  model 
have  been  carried  out  in  studies  of  various  nearshore 
phenomena,  such  as  surf-beat  (Van  Dongeren  et  al., 
1995),  longshore  currents  (Sancho  et  al.,  1995),  infra¬ 
gravity  waves  (Van  Dongeren  et  al.,  1996,  1998;  Van 
Dongeren  and  Svendsen,  2000),  shear  waves  (Sancho 
and  Svendsen,  1998)  and  rip  currents  (Haas  et  al., 
1998;  Svendsen  and  Haas,  1999)  and  the  model  has 
been  compared  to  data  from  the  DELILAH  field 
experiment  (Svendsen  et  al.,  1997). 

The  SHORECIRC  model  was  developed  in  rectan¬ 
gular  Cartesian  coordinates  and  was  only  used  in 
rectangular  domains,  thus  limiting  the  applicability 
of  the  model  to  less  complicated  coastal  environments. 
First,  a  rectangular  grid  is  not  able  to  fit  complicated 
shoreline  boundaries  very  well.  Second,  complicated 
geometries,  such  as  harbors  and  tidal  inlets,  make  the 
uniform-resolution  model  very  expensive  when  a  fine 
grid  is  used.  Therefore,  the  development  of  a  curvilin¬ 
ear  version  of  SHORECIRC  is  necessary  for  its  use  in 
irregular  shaped  domains. 

There  are  numerous  examples  of  curvilinear  grid 
methods  in  the  study  of  waves  and  currents.  Usually 
structured  grid  methods  with  finite-difference  discre¬ 
tizations,  or  unstructured  grid  methods  with  a  finite- 
element  or  finite-volume  approach  are  used  in  model 
developments.  Unstructured  grids  are  more  flexible 
than  structured  grids  to  fit  complicated  boundaries  and 
are  able  to  deal  with  very  complex  geometries.  How¬ 
ever,  the  finite-difference  methods  with  structured 
curvilinear  grids  are  much  simpler  to  program  than 
finite-element  or  finite-volume  methods  and  thus  are 
widely  used  in  fluid  dynamic  fields.  For  example, 
Blumberg  and  Mellor  (1987)  developed  a  3D  coastal 


ocean  circulation  model  (POM)  in  orthogonal  curvi¬ 
linear  coordinates.  Non-orthogonal  boundary-fitted 
grid  models  with  generalized  coordinate  transforma¬ 
tion  were  developed  by  many  authors  (e.g.,  Sheng, 
1986;  Shi  and  Sun,  1995;  Shi  et  al.,  1997)  for  model¬ 
ing  coastal  and  estuarine  processes. 

For  a  generalized  coordinate  transformation,  sever¬ 
al  advantages  of  using  the  contravariant  technique 
have  been  recognized  in  the  derivations  of  hyperbol¬ 
ic-type  equations,  as  shown  by  Sheng  (1986)  and  Shi 
and  Sun  (1995),  among  others.  In  generalized  curvi¬ 
linear  coordinates,  contravariant  and  covariant  compo¬ 
nents  are  two  kinds  of  vector  components  based  on, 
respectively,  a  basis  which  is  locally  tangent  to  the 
curvilinear  coordinate  and  a  reciprocal  basis  which  is 
locally  normal  to  the  curvilinear  coordinate.  The 
designation  “contravariant”  or  “covarianf  ’  technique 
represents  the  choice  of  components  which  are  adop¬ 
ted  as  primitive  variables  in  the  equations  transformed 
from  rectangular  Cartesian  coordinates.  It  has  been 
found  that  both  the  contravariant  technique  and  the 
covariant  technique  are  able  to  simplify  the  trans¬ 
formed  equations  (see,  for  example,  Warsi,  1998),  in 
comparison  to  the  Cartesian  component  method  (see, 
for  example,  Hauser  et  al.,  1985,  1986;  Raghunath  et 
al.,  1987;  Borthwick  and  Barber,  1992).  The  contra¬ 
variant  technique  can  simplify  slip  lateral  boundary 
conditions  compared  to  the  covariant  technique  and  is 
thus  more  conveniently  used  in  hydrodynamic  models 
with  slip  boundary  conditions. 

In  the  present  paper,  a  curvilinear  nearshore  circu¬ 
lation  model  is  developed  based  on  the  quasi-3D 
circulation  equations.  A  generalized  coordinate  trans¬ 
formation  and  the  contravariant  technique  are  used  in 
the  model.  Following  the  Cartesian  version  of  the 
circulation  model  SHORECIRC  (Svendsen  et  al., 
2000),  a  fourth-order  Adams -Bashforth- Moulton 
predictor- corrector  scheme  is  employed  in  the  curvi¬ 
linear  model  to  perform  the  time  integration.  Unlike  the 
spatial  discretization  in  the  Cartesian  version,  we  use  a 
staggered  grid  system  in  the  transformed  image  do¬ 
main.  A  fourth-order  scheme  using  standard  five-point 
finite  differencing  is  used  in  the  first-order  derivative 
terms,  while  a  second-order  scheme  is  used  for  the 
higher  order  derivatives  and  coordinate  metrics  for 
spatial  discretizations.  Various  point  types  are  defined 
in  the  model  code  to  recognize  different  boundary 
conditions,  which  allows  the  model  to  be  used  in 
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complicated  domains  such  as  harbors  and  tidal  inlets. 
We  also  test  the  numerical  convergence,  which  is 
important  for  a  curvilinear  model.  The  model  is  finally 
applied  to  four  cases  involving  either  a  nonorthogonal 
Cartesian  grid  or  generalized  curvilinear  grid.  The  first 
case  is  the  simulation  of  Gourlay’s  (1974)  laboratory 
experiment  in  which  currents  were  generated  on  a 
curved  beach  by  diffracted  waves  in  the  lee  of  a 
breakwater.  We  generated  a  curvilinear  grid  to  fit  the 
curved  beach.  The  breakwater  is  modeled  with  the 
recognition  of  point  types.  The  second  case  is  the 
simulation  of  longshore  circulation  around  a  conical 
island.  A  circular  domain  with  curvilinear  grid  is 
employed  for  this  case.  The  third  case  is  the  longshore 
current  simulations  for  the  longshore  current  experi¬ 
ment  conducted  at  US  Army  Engineer  Research  and 
Development  Center  (Hamilton  and  Ebersole,  2001).  A 
non-rectangular  grid  is  employed  to  fit  the  oblique 
waveguides  used  in  the  physical  experiment.  The  last 
case  is  a  simulation  of  nearshore  circulation  at  Grays 
Harbor,  Washington.  A  curvilinear  grid  is  generated  to 
deal  with  the  complicated  geometry  of  jetties  and  an 
inlet. 

This  paper  is  organized  as  follows.  First,  we  re¬ 
derive  the  circulation  equations  (PS99)  in  generalized 
curvilinear  coordinates.  Next,  model  implementations 
including  numerical  schemes,  vertical  profile  calcula¬ 
tions,  as  well  as  boundary  conditions  are  described. 
Then,  the  numerical  test  on  convergence  and  four 
model  applications  are  carried  out.  The  conclusions 
are  drawn  in  the  final  section. 

2.  Derivations  of  equations  in  curvilinear 
coordinates 

2.7.  Coordinate  transformation 

A  coordinate  transformation  is  introduced  in  the 
general  form 

"  £l  =  ti(xi,x2) 

<  &  =  £,2{xUX2)  (1) 

z  =  z 

\ 

where  (xi,  x2,  z)  are  spatial  independent  variables  in 
rectangular  Cartesian  coordinates  and  (£i,  £2,  z)  are 


new  independent  variables  in  the  transformed  image 
domain,  z  represents  the  vertical  coordinate.  Any 
vector  v  can  be  written  in  Cartesian  coordinates  as 

v  =  vii  +  v2j  T  v3k  (2) 

where  (i,  j,  k)  are  Cartesian  basis  vectors.  In  the 
generalized  curvilinear  coordinates,  the  coordinate 
basis  is  the  combination  of  the  generalized  horizontal 
basis  (a1?  a2)  and  the  vertical  Cartesian  base  vector  k 
and  is  written  as  (a1?  a2,  k).  Using  the  new  basis  (ai, 
a2,  k),  the  vector  v  can  be  described  as 

v  =  vaaa  +  v3k  (a  =1,2)  (3) 

where  va  are  contravariant  horizontal  components  of 
the  vector;  v3  is  the  vertical  Cartesian  component  of  the 
vector. 

The  relationship  between  horizontal  components  in 
the  Cartesian  and  contravariant  components  in  the 
curvilinear  coordinates  may  be  written  as 

V*  =  |^v/,  («,/?=  1,2)  (4) 

oxp 

where  vp  represents  the  Cartesian  components.Follow- 
ing  the  velocity  splitting  method  in  PS99,  the  contra¬ 
variant  components  of  instantaneous  horizontal 
velocity  can  be  written  as 

=  u'«  +  <,  +  Va  +  VI  (5) 

where  i/a,  w“,  Ea,  and  V\  are,  respectively,  the  turbu¬ 
lence  component,  the  wave  component,  the  component 
of  depth-averaged  and  short-wave-averaged  velocity, 
and  the  vertical  variation  of  the  short-wave-averaged 
velocity. 

2.2.  Depth-integrated,  short-wave-averaged  equations 

The  model  equations  in  PS99  were  derived  in 
rectangular  Cartesian  coordinates  and  described  in 
terms  of  Cartesian  tensor  representations.  It  may  not 
be  appropriate  if  we  only  transform  the  final  Eqs.  (4) 
and  (45)  in  PS99  from  Cartesian  coordinates  into 
curvilinear  coordinates  since  use  of  the  direct  trans¬ 
formation  does  not  allow  a  simplification  of  the 
equations  using  the  contravariant  technique.  In  addi¬ 
tion,  we  can  not  guarantee  that  the  curvilinear  counter- 


102 


F.  Shi  et  al.  /  Coastal  Engineering  49  (2003)  99-124 


parts  of  the  3D  coefficients  M,  A,  D  and  B  in  PS99)  are 
in  correct  form  after  the  generalized  coordinate  trans¬ 
formation.  Therefore,  we  re-derive  the  equations  in 
generalized  coordinates,  starting  with  the  depth-inte¬ 
grated,  short-wave-averaged  equations  as  shown  in 
Eqs.  (4)  and  (5)  in  PS99. 

The  depth-integrated,  short-wave-averaged  conti¬ 
nuity  equation  (PS99,  Eq.  (4))  in  terms  of  contravariant 
components  is  given  by 


ac 

a  t 


l 


He 


W&V*h)  =  0 


(6) 


where  g0  is  the  determinant  of  the  metric  tensor  gap, 


gn 

gl2 

g2l 

g22 

in  which 


(7) 


t  s  and  Tb  are  the  contravariant  components  of  the 
surface  shear  stress  and  the  bottom  shear  stress, 
respectively.  S and  Tx(j  are  the  contravariant  compo¬ 
nents  of  the  short-wave-induced  radiation  stress  tensor 
and  the  depth-integrated  Reynolds’  stress  tensor,  re¬ 
spectively.  ()$  presents  the  covariant  derivative  of  a 
second-order  tensor,  given  by 


3 

Up 


_|_  rpyPp 


a 

yP 


(11) 


where  rfp  is  Christoffel  symbol  of  the  second  kind. 
The  detailed  derivation  of  Eq.  (9)  and  calculations  of 
Tffi  and  S°jf_  can  be  found  in  Appendix  A. 

'  In  (9),  ’FhV?V?<lz+  Q^vfiO  +  VfiOQd  repre- 
sents  the  effects  of  the  vertical  nonuniformity  of  the 
short-wave-averaged  velocities  and  give  rise  to  the 
dispersive  mixing  as  described  in  PS99.  These  terms 
are  evaluated  using  the  solution  for  Vf. 


0  Xy  3  Xy 

8x11  =  JC.Wp 


(8) 


To  get  Eq.  (6),  we  start  with  the  general  form  of  the 
continuity  equation  and  then  write  it  in  a  tensor- 
invariant  form  that  can  be  simplified  to  the  contra¬ 
variant  component  equation.  The  detailed  derivation 
of  Eq.  (6)  can  be  found  in  Appendix  A. 

Similarly,  the  depth-integrated,  short-wave- aver¬ 
aged  momentum  equation  can  be  written  in  terms  of 
contravariant  components  as 


2.3.  Solution  for  Vf 

To  get  the  solution  for  Vf,  we  start  with  the 
horizontal  momentum  equation  as  in  PS99.  The  fol¬ 
lowing  tensor-invariant  form  of  the  equation  is  ob¬ 
tained  from  the  general  form  of  the  Eq.  (55)  shown  in 
Appendix  A. 


YtiFh) 


nmi 


vaPh  +  [  vfvfdz  +  Qrv?(t) 

J-ho 

1 


T  UP 

CP  ' 


1  ,,fi  ,  3{ 

7s'  +skW/ 


.Pot 


--<+-^b=0 

P  P 


(9) 


where  Qf  represents  the  contravariant  components  of 
short-wave-induced  volume  flux  as  defined  in  PS99. 
are  the  contravariant  metric  defined  as 


where  p  is  the  instantaneous  pressure;  w  is  the  vertical 
component  of  the  instantaneous  velocity.  Introducing 
Eq.  (5)  into  Eq.  (12)  and  averaging  over  a  wave  period 
leads  to 

2E  +  +  ( y«vf  +  vaP  +  v«P  +  v[va- 

dt  dt  v  1  1  1  1 

+  ui  +  IFiF]  +  ~-(v?iv  +  vxw 

/  oz  \ 

+  +  =  -L^E-gP*  (13) 

7  P  OQp 


0Xy  0Xy 


(10) 


where  W represents  the  vertical  component  of  the  time 
averaging  velocity.  Assuming  hydrostatic  pressure  and 
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the  eddy  viscosity  closure  (see  Appendix  A),  we  get 
the  tensor-invariant  form  of  PS99  (Eq.  (16)) 


dt  dz 


dF[ 

dz 


-  ( v*vh>  +  viFp  +  v\vln  +  w 


a  y* 

dz 


+  [MgyPvxy+gyxv%  +  Yz(Vtgl>xjf) (14) 


where 


By  comparing  Eq.  (17)  against  PS99,  Eq.  (19),  the 
major  changes  of  the  equation  after  the  transformation 
can  be  found  that  the  regular  horizontal  derivatives  in 
the  Cartesian  equation  become  the  covariant  deriva¬ 
tives  in  the  generalized  curvilinear  equation.  The 
vertical  derivative  terms  in  the  generalized  equation 
appear  to  be  in  the  same  form  as  that  in  the  Cartesian 
equation  since  there  is  no  transformation  in  the  vertical 
direction.  Therefore,  the  solution  for  Vf  should  be  in 
the  same  form  as  in  Cartesian  coordinates  as  described 
below. 

After  the  non-dimensional  analysis  and  perturba¬ 
tion  as  in  PS99,  we  use  a  perturbation  expansion  of 
the  type 


and  Q/j  represents  the  covariant  derivative  of  a  first- 
order  tensor  defined  by  Eq.  (43). Using  the  continuity 
equation  (Eq.  (6)),  the  depth-integrated  momentum 
Eq.  (9)  may  be  written  as 


dt 


1 


1 


1 


~  ~  _a  _  ~  rjiV.fi _ ■ 

ph  ’P  ph  B  ph  ’P  h 


x  V?V?dz  +  Ql vf(0  +  V?(0Qij  (16) 


Eq.  (16)  can  be  used  to  rewrite  Eq.  (14)  as 

f)  ytx.  a 

^--(vA) 

dt  dzy  dz  J 
_  (  1  oajS  _  fa  |  TB 

-{ph’*  1  +  ph 


-  ( +  vfvfp  +  vfvfp  +  w^l 


i  (  F 


dz 


V?vfdz  +  QlvfCO  +  V?mi 


h  \J-h 

{ [vt(gypv«+gy*v% 


+yAv‘8  WJ  Fh  F 


V?  =  v(0)  +5F1a(1)  +  ...  (18) 

where  6  (  ~  0.1)  represents  the  size  of  the  short- wave- 
induced  quantities;  Vi(0)  and  fT(l)  represent  zero- 
and  first-order  of  the  contravariant  components.  We 
only  consider  the  first  two  terms  on  the  right  hand  side 
of  (18)  in  the  present  paper.  The  governing  equation 
for  Fia(0)  is 


9<0)  3  /  3F“(0)\  _F< 

3 1  dz  y  dz  J 

where 


1 

ph 


Fa  =— ST  -/a  +  ^- 

,p  ph 


(19) 


(20) 


The  solution  for  Eia(0)  will  be  the  same  as  the 
solution  in  Cartesian  coordinates  (PS99,  Eqs.  (34)  and 
(37))  except  the  solution  value  represents  the  contra¬ 
variant  component  of  the  vector. 

Similarly,  the  solution  for  Eja(0)  can  be  written  in 
the  same  form  as  in  PS99,  Eq.  (42)  except  that 


(z,  t)  =  -  (  +  vf{0)V«p  +  Kf(0)F“J 


(o) 

ft 


dV\ 
+  W  1 


a(0)  ' 


dz 


(17) 


+  Qlv(°\o  +  v(°\Q  Qi 


(21) 
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2.4.  Results  for  the  integrals 

The  derivation  for  the  integrals  required  in  Eq.  (9)  is 
very  similar  to  that  in  PS99,  noticing  that  all  the 
horizontal  derivatives  in  PS99  should  be  replaced  by 
the  corresponding  covariant  derivatives  in  generalized 
coordinates. 

In  PS99,  Eq.  (94)  was  used  to  replace  W  in  PS99, 
Eq.  (93).  The  representation  of  W  in  general  coordi¬ 
nates  may  be  written  as 


The  final  tensor  form  of  the  momentum  Eq.  (16) 
(similar  to  PS99,  Eq.  (45))  can  be  written  as 

l~t  ( Vxh)  +  ( F  Ph  +  Aaps F)p  +  3  (Sali  +  PMap)j 

+  ghgPa -S  +  f+  ~  h(DslsV«  +  DdxV% 
oQp  p 

-  =  0  (25) 


w  =  - 


{F  +  vC\-k))2f 

dxp 


(ho  +  z)V^  +  f 

J-ho 


Vmdz 
v \$  uz 


where  the  tensors  A,  B,  D  and  M  have  the  same 
definitions  as  in  PS99  except  all  the  horizontal  deriv¬ 
atives  of  should  be  replaced  by  the  corresponding 
co variant  derivatives. 

Expanding  the  covariant  derivatives  in  Eq.  (25), 
the  momentum  equations  in  generalized  curvilinear 
coordinates  can  also  be  written  as 


In  the  derivation  of  Eq.  (22),  we  use  the  following 
relations: 


JA±+V>  <-*„)!! 


rr^dz 


(23) 


2(^)  +  _ L3L 


[^{PPh+AaltdP)\ 


+  ( P  Ph + AyPS  P)r;p  +  -2—  A. 

yP  PVTq  Up 

X  [v®*-''  +  +  -(«"’  +  pM,t)r. 


yp 


-  h(DSflFs+DSoiPs)]}  +  [TyP  -  h(DspV} 

+ DSyPs)]r;p  ~J=  A.  [^(hB'pvft] 

-  0 hBypvfrr^,  =  o 


(26) 


Eqs.  (6)  and  (26)  are  the  nearshore  circulation 
equations  in  generalized  curvilinear  coordinates. 


Following  the  steps  in  PS99  leads  to  the  results  for 
the  integrals  in  Eq.  (16)  in  general  coordinates: 

f  pvfdz  +  Q*wVf( 0  +  V^OQt 

J  ~ho 

=  Maf}  +  Aaed  Vs  -  h(Dsp Fs  +  DSx Ps  +  BaP  Ps ) 

(24) 


3.  Model  implementation 

3.1.  Numerical  scheme 

Following  the  Cartesian  version  of  the  SHORE- 
CIRC  model,  the  governing  equations  (Eqs.  (6)  and 
(26))  are  solved  using  the  fourth-order  Adams -Bash- 
forth-Moulton  predictor- corrector  scheme  to  per- 
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form  the  time  integration.  A  fourth-order  scheme 
using  standard  five-point  finite  differencing  is  used 
for  the  first-order  spatial  derivative  terms  and  a 
second-order  central  scheme  is  used  for  higher  order 
derivatives,  metric  tensors  and  Christoffel  symbols.  In 
contrast  to  the  non-staggered  grid  used  in  the  Carte¬ 
sian  version,  we  employ  a  staggered  grid  in  the  (£i, 
£2)  plane  with  which  the  numerical  noise  level  was 
found  to  become  lower  than  that  with  a  non-staggered 
grid.  The  staggered  grid  arrangement  and  the  detailed 
numerical  schemes  can  be  found  in  Appendix  C.  In 
addition,  various  point  types  are  defined  on  the 
staggered  grid  to  recognize  different  boundary  con¬ 
ditions  and  to  deal  with  complicated  structures  in  a 
computational  domain.  No  numerical  filtering  is  need¬ 
ed  in  the  code. 

3.2.  Calculation  of  vertical  current  profiles 

As  described  in  Eq.  (18),  the  depth- varying  current 
velocity  is  split  into  two  parts.  The  first  part  V i(0)  is 
primarily  the  component  generated  by  the  local 
external  forcing  F*  in  Eq.  (20)  while  the  second  part 
n(1)  is  generated  by  the  advective  terms  shown  in 
Eq.  (17).  After  a  scale  analysis  in  PS99,  a  reasonable 
assumption  is  made  that  the  first  component  is  much 
larger  than  the  second  one.  In  addition,  the  expres¬ 
sions  of  coefficients  A,  B ,  D  and  M  in  PS99  also  show 
that  the  contributions  to  these  coefficients  from 
can  be  expressed  in  terms  of  the  f/0)  component  in 
the  values  of  the  current- current  and  current- wave 
interaction  terms.  Therefore,  we  only  use  to 

present  the  vertical  variation  of  Vf.  E/0)  is  given  by 


—  d\^2  +  +/la  3-fla 

(27) 

where 

Z=z  +  h 

(28) 

in  which  h  is  the  water  depth  from  the  mean  water 
level  to  the  bottom.  In  Eq.  (27) 


streaming  and  fa  represents  the  short-wave  forcing 
(see  Putrevu  and  Svendsen,  1995  for  detail). 
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Substituting  Eq.  (27)  into  the  integration  forms  of  the 
coefficients  A,  B,  D  and  M  in  PS99,  we  can  get  the 
final  expressions  of  the  dispersive  mixing  coefficients 
as  shown  in  Appendix  B. 


3.3.  Boundary  conditions 


There  are  several  types  of  boundary  conditions 
implemented  in  the  curvilinear  model.  First,  slip 
boundary  conditions  can  be  used  for  lateral  walls, 
shorelines  and  boundaries  of  structures  inside  a  com¬ 
putational  domain.  As  mentioned  in  the  introduction  of 
the  paper,  the  contravariant  technique  can  simplify  the 
expressions  of  slip  boundary  conditions  in  curvilinear 
coordinates.  For  a  curvilinear  boundary,  the  slip 
boundary  condition  can  be  simply  expressed  as 

Fa  =  0.  (32) 

Second,  a  specified-flux  boundary  condition  is 
implemented  by  using  the  contravariant  component. 
The  contravariant  component  of  the  specified  velocity 
(or  flux)  can  be  obtained  using  the  transformation 
relation  (Eq.  (4)). 

Third,  a  periodic  boundary  condition  along  cross¬ 
shore  boundaries  at  the  two  ends  of  the  domain  is  also 
an  option  in  this  model  for  simulations  of  uniform 
beaches.  The  implementation  of  the  periodic  bound¬ 
ary  condition  is  similar  to  that  in  the  Cartesian  version 
of  the  SHORECIRC  model  except  that  all  vector 
variables  are  taken  in  contravariant  forms.  The  peri¬ 
odic  boundary  condition  requires  both  the  bathymetry 
and  grid  to  be  periodic. 


j_((v-sr+T“-T“)-/, 


(29) 


4.  Numerical  test  and  application 

4.1.  Convergence  test 


where  ( V  -S)a  is  contravariant  component  of  V  -S;  ts*  Convergence  is  a  very  important  numerical  proper- 

represents  the  shear  stress  associated  with  the  steady  ty  for  a  curvilinear  model  since  variable  grid  spacing  is 
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generally  used  in  a  grid  system.  As  a  simple  test  case 
(Shi  et  al.,  2001),  the  evolution  of  waves  in  a  rectan¬ 
gular  basin  is  calculated  to  test  the  convergence  with 
both  space  and  time  discretization. 

The  basin  dimensions  are  20  x  20  m,  and  the  water 
depth  is  0.5  m  constant  over  the  basin.  The  initial 
condition  is  provided  by  a  motionless  Gaussian  hump 
of  water  with  its  center  located  at  the  center  of  the 
basin  (xc,  yc): 

C(x,y,t=0)=H0exp{-y[(x  -  xc)2  +  (y  -  ycf}/L2}, 

(33) 


elevations  at  t=  10  s.  The  RMS  difference  of  surface 
elevation  is  defined  by 


RMS/;  = 


N 


-  £p+ i(v))2 

j=  1  7=1 

m  x  n 


(36) 


where  CP  represents  the  calculated  surface  elevations 
in  the  case  p.  ( m ,  n)  are  the  grid  dimensions  in  the 
case  of/?  =  8.  Fig.  1  shows  that  the  logarithmic  RMS 
differences  decrease  linearly  when  grid  spacing  de¬ 
crease.  Then,  the  Cauchy  convergence  rate  defined 
by  the  following  formula  can  be  calculated. 


R 

T? 

II 

o 

II 

(34) 

v(x,y,t  =  0)  =  0, 

(35) 

where  H0  is  the  initial  height  of  the  hump;  L  is  a  scale 
length;  y  is  the  shape  coefficient,  and  (pcc,  yc)  is  the 
coordinate  at  the  center  of  the  domain.  We  chose 
//o  =  0.2  m,  L=  1  m,  y  =  0.4,  and xc=yc  =  10  m. 

To  test  the  grid  convergence,  we  keep  the  time 
step  as  a  constant,  i.e.,  A^=0.01  s,  and  adopt  a 
sequence  of  different  grid  spacing  Ax  x  /?,  where 
Ax  =  0.05  m  and  p=  1,2,. .  .,8  (The  Courant  Numbers 
are  less  than  0.5  in  all  the  cases).  Fig.  1  shows  the 
convergence  rate  with  grid  refinement  that  is  demon¬ 
strated  by  the  RMS  differences  of  simulated  surface 


Fig.  1.  Convergence  rates  with  grid  refinement. 


log(RMSJ,/RMSj>+1) 

\og(Axp/Axp+1)  [  1 

The  averaged  R  calculated  in  this  case  is  3.58  which 
is  very  reasonable  since  both  fourth-order  and  sec¬ 
ond-order  schemes  are  utilized  in  the  model. 

Similarly,  the  convergence  with  time  discretization 
is  examined  by  using  a  sequence  of  time  steps  from 
0.004  to  0.02  s  and  keeping  a  constant  grid  spacing 
Ax  =  0.2  m.  The  convergence  rate  with  time  step 
refinement  is  shown  in  Fig.  2.  The  averaged  conver¬ 
gence  rate  is  2.32  which  is  a  little  lower  than  the  grid 
spacing  convergence  rate.  It  is  found  that  the  time 
refinement  convergence  rate  can  be  improved  to 
about  3.0  by  using  an  under-relaxed  iteration  tech¬ 
nique  (Wei  and  Kirby,  1995)  during  the  corrector 
stage  in  the  code.  We  do  not  use  the  iteration 
technique  for  computational  economy  and  assume 
the  convergence  rate  without  iteration  to  be  sufficient 
for  the  model. 

4.2.  Simulation  of  laboratory  experiment  of  Gourlay 

Gourlay  (1974)  carried  out  a  laboratory  experi¬ 
ment  in  which  currents  were  generated  on  a  curved 
beach  by  diffracted  waves  in  the  lee  of  a  breakwater. 
The  purpose  of  the  experiment  was  to  demonstrate 
that  longshore  variation  of  the  wave  height  could 
generate  longshore  currents.  Fig.  3  shows  the  labo¬ 
ratory  setup  for  Gourlay’s  experiment.  A  1  on  10 
concrete  beach  is  parallel  to  the  incoming  wave  crests 
in  the  exposed  zone.  In  the  shadow  zone,  the  slope  is 
also  1/10  towards  the  curved  beach  which  has  a 
constant  radius  centered  on  the  breakwater  tip.  The 
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wave  basin  was  thus  designed  so  that  the  shoreline 
was  everywhere  approximately  parallel  with  the  dif¬ 
fracted  wave  crests. 


Fig.  3.  Layout  of  Gourlay’s  experiment. 


Winer  (1988)  employed  a  numerical  model  to 
simulate  Gourlay’s  laboratory  experiment.  He  used 
REF/DIF-1  (Kirby  et  al.,  2002)  as  the  wave  driver 
and  used  depth-averaged  shallow  water  equations 
which  include  radiation  stress  terms  in  the  circulation 
model.  Because  a  Cartesian  coordinate  system  was 
adopted  in  the  numerical  model,  a  stair-type  boundary 
was  employed  at  the  curved  beach.  This  boundary 
effect  was  shown  in  the  contour  plots  of  both  wave  set¬ 
up  and  velocity  amplitude  in  his  results. 

We  use  the  present  curvilinear  model  to  simulate  the 
experiment.  A  curvilinear  grid  shown  in  Fig.  4  is  used 
to  fit  the  curved  shoreline  boundary.  Rather  than  the 
small  computational  domain  Winer  chose  from  the 
shoreline  to  3  m  offshore  of  the  breakwater,  a  larger 
computational  domain  including  the  wave  paddle  area 
is  employed  to  provide  enough  reservoir  to  feed  the 
wave  set-up,  as  in  the  laboratory  experiment.  The  grid 
sizes  are  about  0.1  m  in  the  offshore  region  and  the 
exposed  straight  shoreline  region  and  less  than  0.1  m 
around  the  curved  beach. 

REF/DIF-1  is  used  as  the  wave-driver  in  a  subrou¬ 
tine  of  the  model  code  and  provides  the  circulation 
model  with  radiation  stresses,  short- wave-induced 
volume  flux,  breaking  wave  energy  dissipation,  and 
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Fig.  5.  Wave  set-up  contours  from  Gourlay’s  experiment  (from  Gourlay,  1974). 


wave  bottom  velocity.  Because  REF/DIF-1  is  operat¬ 
ing  in  a  Cartesian  grid,  an  interpolation  is  needed  for 
the  data  transfer  from  the  Cartesian  grid  to  the  curvi¬ 
linear  grid  (see  Appendix  D  for  detail).  In  REF/DIF-1, 
the  wet  and  dry  grid  technique  proposed  by  Kirby  and 
Dalrymple  (1986)  is  used  for  treating  the  dry  grid 
points  as  though  they  have  a  very  small  depth  of  water 
(1  mm  in  the  present  paper).  A  wave  height  of  9.1  cm 
with  a  wave  period  of  1 .5  s  is  used  in  REF/DIF- 1  as  the 
incident  wave  conditions,  as  in  Winer’s  simulation. 

The  bottom  stress  formulation  (Svendsen  et  al., 
2000)  with  a  constant  bottom  friction  coefficient 


fcw  =  0.008  is  adopted  in  the  case.  The  turbulence 
mixing  coefficients  in  the  eddy  viscosity  formulation 
are  chosen  as  the  same  values  as  in  Svendsen  et  al. 
(2000).  The  coefficients  are  also  used  in  the  following 
laboratory  experiment  cases. 

Figs.  5  and  6  show  the  experimental  results  for  the 
mean  water  surface  contours  and  the  corresponding 
results  obtained  from  the  numerical  model,  respective¬ 
ly.  The  comparison  demonstrates  that  the  numerical 
results  agree  very  well  with  the  experiment  results.  In 
contrast  to  Winer’s  results,  which  failed  to  show  any 
set-up  in  the  area  of  the  connection  of  the  shoreline  and 


Fig.  6.  Wave  set-up  contours  from  numerical  results. 
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Fig.  7.  Contours  of  current  velocities  and  streamlines  from  Gourlay’s  experiment  (from  Gourlay,  1974). 


the  breakwater,  the  present  model  results  show  a  wave 
set-up  in  this  region  which  is  consistent  with  the 
laboratory  results. 

To  make  model/data  comparisons  for  the  circula¬ 
tion,  we  present  Gourlay’s  experimental  results  in  Fig. 
7.  The  figure  illustrates  both  the  contours  for  the 
magnitude  of  the  velocity  and  some  streamlines  indi¬ 
cating  the  direction  of  the  flow.  The  corresponding 
results  from  the  numerical  model  are  shown  in  Fig.  8. 
It  can  be  seen  from  the  comparison  that  the  overall 
agreement  is  fairly  good.  Both  the  magnitude  and  the 
locations  of  the  maximum  currents  obtained  in  the 


numerical  model  are  similar  to  that  in  the  laboratory 
experiment.  The  shoreward  oriented  current  in  the  area 
directly  exposed  to  the  incident  waves  is  very  well 
shown  in  the  numerical  results.  The  present  model  also 
predicts  the  location  of  the  primary  eddy  center  and 
current  magnitude  behind  the  breakwater  as  shown  in 
Fig.  8.  However,  the  model  fails  to  predict  the  stagna¬ 
tion  eddy  in  the  comer  behind  the  breakwater  although 
it  does  provide  the  set-up  variation  in  that  area  as 
shown  in  Fig.  6.  The  parabolic  wave  model  is  sus¬ 
pected  responsible  to  cause  the  discrepancy  because  it 
is  not  appropriate  model  in  predicting  wave  diffraction 


Fig.  8.  Contours  of  current  velocities  and  velocity  vectors  from  numerical  results. 
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periodic  boundaries 


Fig.  9.  Computational  grid  for  simulation  of  longshore  circulation 
around  a  conical  island. 

in  the  sheltered  region.  Further  studies  using  other 
wave  models  may  be  needed  to  precisely  simulate  this 
phenomena. 

4.3.  Longshore  circulation  around  a  Conical  Island 

Mei  and  Angelides  (1977)  conducted  a  semi-theo¬ 
retical  study  of  the  averaged  circulation  caused  by 
waves  breaking  on  the  constant  slope  beach  around  a 
circular  island  in  a  sea  of  constant  depth.  They  used  the 


geometric  refraction  method  and  the  energy  conserva¬ 
tion  theory  for  wave  calculations.  For  calculations  of 
circulation,  they  used  very  simple  depth-and  time- 
averaged  equations  of  motion  that  include  pressure 
gradient,  radiation  stress  and  bottom  friction  terms. 
Polar  coordinates  were  used  in  the  circulation  model. 
They  concluded  that  the  wave  refraction  pattern  deter¬ 
mines  the  extent  of  longshore  circulation.  In  particular, 
for  an  island  so  large  that  there  is  a  lee  shore  which  is 
not  affected  by  the  refracting  waves,  two  current  cells 
are  found  on  two  sides  of  the  island  with  respect  to  the 
wave  direction,  without  significant  penetration  in  the 
lee  shore  region.  This  conclusion  was  used  to  explain 
the  creation  of  two  sand  spits  for  abrasive  island  shores. 

As  a  simple  and  efficient  case  for  testing  the 
coordinate  transformation  within  a  circular  domain, 
the  so-called  Targe  island’  case  is  carried  out  in  the 
present  paper  as  in  Mei  and  Angelides  (1977).  The 
radius  of  the  island  is  3048  m  (10,000  ft).  The  water 
depth  of  the  flat  seabed  is  30.48  m  (100  ft)  and  the 
beach  slope  we  adopted  in  the  simulation  is  1/40.  The 
incident  wave  height  and  wave  period  are,  respective¬ 
ly,  6  m  and  10  s.  REF/DIF-1  is  again  used  for  wave 
calculation  and  the  wet  and  dry  grid  technique  is  also 
used  to  deal  with  the  island  geometry. 

Fig.  9  shows  the  computational  grid  generated 
according  to  the  following  formula: 

x(i,j)  =  [3048  +  15 (/  -  1)  +  0.865(z  -  l)2]cos(/) 

_y(/j)  =  [3048+15(/—  l)+0.865(z-  l)2]sin(/)  (38) 

where  x(i,  j )  and  y(i,  j )  are  grid  point  coordinates  at 
(/,/);  i=  1,2,. .  .,52  and  j=  1,2,. .  .,361  (where  the  argu- 


Fig.  10.  Calculated  wave  set-up. 
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Fig.  11.  Vectors  of  volume  flux  calculated  from  the  numerical  model. 


ment  of  the  sine  and  cosine  is  assumed  to  be  degrees 
in  Eq.  (38)).  The  figure  shows  a  finer  grid  used  in  the 
nearshore  region  to  resolve  the  surf  zone.  The  mini¬ 
mum  grid  spacing  is  15.8  m  in  the  nearshore  region 
and  the  maximum  is  60  m  in  the  offshore  region. 


Notice  that  there  are  four  boundaries  in  the  computa¬ 
tional  domain,  rather  than  two  boundaries  usually 
used  in  the  polar  coordinates.  Besides  the  island 
boundary  and  the  offshore  boundary  shown  in  Fig. 
9,  there  are  two  lateral  boundaries  linked  on  the  lee 


Fig.  12.  Plan  view  of  the  LSTF. 
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Fig.  13.  Computational  grid. 


side  along  a  island-normal  coordinate  line.  Periodic 
boundary  conditions  are  used  at  the  two  lateral 
boundaries. 

Fig.  10  shows  the  mean  surface  elevation  calcu¬ 
lated  from  the  circulation  model.  It  is  clear  that  a 
wave  set-up  occurs  at  the  front  shore  and  a  slight 
wave  set-down  in  the  shoaling  zone.  The  current 
pattern  can  be  demonstrated  by  the  volume  flux  of 
the  calculated  circulation  shown  in  Fig.  11.  There  are 
strong  longshore  currents  along  the  front  shore.  The 
currents  turn  back  around  the  edges  of  the  lee  shore 
and  complete  the  two  cells  through  the  shoaling  zone. 
This  result  is  consistent  with  the  results  of  Mei  and 
Angelides  (1977).  The  phenomenon  was  thus  used  to 


explain  the  formation  of  spits  protruding  from  the 
island. 

4.4.  Longshore  current  simulation  in  an  obliquely 
quadrilateral  domain 

Recently,  a  longshore  current  experiment  was 
carried  out  in  the  Large-scale  Sediment  Transport 
Facility  (LSTF)  at  the  U.S.  Army  Engineer  Research 
and  Development  Center’s  Coastal  and  Hydraulics 
Laboratory  (Hamilton  and  Ebersole,  2001).  The  LSTF 
has  a  concrete  beach  which  has  a  longshore  dimen¬ 
sion  of  31  m  and  a  cross-shore  dimension  of  21  m, 
with  a  plane  slope  of  1:30.  To  reproduce  a  longshore 


Fig.  14.  Depth-averaged  current  velocity  field. 
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uniform  current  in  the  finite-length  wave  basin,  20 
independent  pumps  and  channels  are  used  to  control 
the  cross-shore  distribution  of  the  mean  longshore 
current.  To  minimize  wave  diffraction  into  the  flow 
channels,  two  obliquely  oriented  waveguides  oriented 
10°  from  the  shoreline-normal  direction  were  set  on 
lateral  boundaries.  Fig.  12  shows  a  plan  view  of  the 
LSTF  and  locations  of  waveguides  and  measurement 
transects  (dashed  lines).  In  the  figure,  Qp  is  the  total 


longshore  flow  rate  actively  pumped  through  the 
external  recirculation  system.  If  the  wave-induced 
longshore  current  does  not  match  Qp ,  an  internal 
recirculation  Qr  will  develop  as  shown  in  the  figure. 
The  Qs  is  the  total  longshore  flow  rate  between  the 
wave  set-up  limit  and  the  point  of  transition  where  the 
mean  longshore  current  reverses  direction.  In  the 
physical  experiment,  the  discharge  from  the  pumps 
was  adjusted  to  minimize  the  internal  recirculation  Qr. 


Fig.  15.  Longhore  current  comparisons  at  the  measurement  transects  (circles:  measurement,  solid  lines:  numerical  results). 
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This  technique  was  previously  used  by  Visser  (1980, 
1984,  1991). 

For  the  simulation  of  a  perfect  longshore  uniform 
current  without  internal  recirculation  Qr ,  the  Cartesian 
version  of  SHORECIRC  may  be  a  proper  model  for 
longshore  current  simulations  because  there  is  no 
boundary  effect  caused  by  the  oblique  waveguides. 
For  a  more  general  case,  however,  the  internal  recir¬ 
culation  exists  and  the  boundary  effect  of  oblique 
waveguides  shows  up  near  the  lateral  boundaries.  A 
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cross-shore  location  x  (m) 


boundary-fitted  grid  model  is  more  appropriate  for  this 
case  than  a  rectangular  grid  model.  A  boundary-fitted 
grid  is  generated  as  shown  in  Fig.  13.  The  grid  sizes 
are  0.5  m  in  longshore  direction  and  0.26  m  along  the 
waveguide  direction  (0.25  m  between  two  adjacent 
longshore-direction  grid  lines). 

Again  we  use  REF/DIF-1  as  the  wave-driver  in  a 
subroutine  of  the  curvilinear  SHORECIRC  model.  As 
a  specified  volume  flux  boundary  condition,  the  mea¬ 
sured  flux  at  the  lateral  boundaries  is  employed  by 
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Fig.  16.  Mean  water  level  comparisons  at  the  measurement  transects  (circles:  measurement,  solid  lines:  numerical  results). 
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interpolating  the  data  into  the  model  boundary  points. 
The  offshore  boundary  and  shoreline  boundary  are 
vertical  walls  in  this  case.  A  regular  wave  case  was 
simulated  from  the  physical  experiment  with  a  wave 
period  of  2.5  s,  0.182  m  incident  wave  height  and  10- 
degree  wave  angle.  The  water  depth  decreases  linearly 
from  0.667  m  at  the  offshore  boundary  to  0.002  m  at 
the  shoreline  boundary.  The  time  step  used  in  this  case 


is  0.05  s.  Wave -current  interaction  is  taken  into 
account  in  the  model  by  calling  the  wave-driver 
subroutine  every  10  s.  The  model  reached  steady  state 
after  200  s. 

Fig.  14  shows  the  calculated  depth-averaged  cur¬ 
rent  field.  It  shows  that  the  steady  state  longshore 
current  appears  approximately  longshore  uniform.  The 
peak  currents  are  around  x  =  12-15  m  in  the  surfzone 
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Fig.  17.  Comparison  of  the  vertical  profile  of  cross-shore  current  at  Y27  (circles:  measurement,  solid  lines:  numerical  results). 
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Fig.  18.  Bathymetry  of  Grays  Harbor  (depths  below  mean  tide  level 
in  meters). 

while  the  breaking  line  is  around  x  =  9  m  in  both  the 
physical  experiment  and  the  REF/DIF-1  results.  There 
is  a  very  weak  internal  recirculation  outside  the  surf- 
zone,  found  in  the  numerical  results  and  confirmed  by 
the  experimental  data. 

The  longshore  current  comparisons  are  made  be¬ 
tween  the  numerical  results  and  experiment  measure¬ 
ments  at  eight  measurement  transects  shown  by  dashed 
lines  in  Fig.  12.  Fig.  15  shows  the  comparisons,  where 
y=  5.7,  9.7,  13.7,  17.7,  21.7,  25.7,  29.7  and  30.7  m  in 
the  figure  are  the  model  coordinates  and  represent  the 
locations  of  the  corresponding  measurement  transects 
Y39,  Y35,  Y31,  Y27,  Y23,  Y19,  Y15,  and  Y14, 
respectively,  shown  in  Fig.  12.  As  shown  in  Fig.  15, 
the  current  amplitudes  and  the  locations  of  peak 
currents  are  well  predicted.  Some  apparent  disagree¬ 
ments  are  found  in  the  downstream  transects  of 
y  =  29.1  and  30.7  m,  where  the  measured  current  data 


show  a  different  cross-shore  variation  around  x  =  12  m. 
This  particular  feature  is  not  resolved  by  the  model 
predictions.  Also,  an  under-prediction  of  the  currents  is 
found  near  the  shoreline  region  in  the  same  two  trans¬ 
ects.  In  Fig.  16,  we  compare  the  measured  mean  water 
levels  and  the  numerical  results.  The  figure  shows  that 
very  good  agreement  is  obtained  from  the  comparisons 
of  wave  set-up  and  set-down. 

The  FSTF  measurements  also  provide  information 
about  the  vertical  variation  of  the  currents.  As  de- 


Fig.  19.  Computational  grid  for  the  Grays  Harbor  case  (Maximum 
grid  size  is  296.4  m  and  minimum  is  59.1  m). 
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scribed  in  the  model  theory,  the  quasi-3D  circulation 
model  is  able  to  predict  the  vertical  profile  of  wave- 
induced  current.  Using  Eq.  (27),  we  can  easily  calculate 
the  vertical  distributions  of  cross-shore  and  alongshore 
velocity.  Fig.  17  shows  the  model/data  comparisons  of 
cross-shore  velocities  at  measurement  points  in  transect 
Y27.  It  is  shown  that  most  of  the  numerical  results 
agree  well  with  the  measurement  data  except  the  results 
at  x  =  10.9  m  where  significant  deviations  occur.  As 
shown  by  Svendsen  and  Putrevu  (1994)  the  vertical 
variations  of  the  current  are  instrumental  in  providing 
the  current  lateral  mixing. 

4.5.  Simulation  of  nearshore  circulation  at  Grays 
Harbor 

Grays  Harbor  is  a  jettied  entrance  on  the  Washington 
coast,  located  approximately  30  km  north  of  Willapa 


Bay.  We  choose  Grays  Harbor  as  the  computational 
domain  in  the  case  study  of  the  curvilinear  version  of 
SHORECIRC  because  it  includes  irregular  coast  lines, 
jetties,  and  complicated  geometry.  The  case  illustrates 
the  versatility  of  the  curvilinear  model  to  deal  with 
complicated  geometries. 

Fig.  18  shows  the  bathymetry  of  Grays  Harbor. 
There  are  two  jetties  located  on  the  south  and  north 
side  of  the  harbor  entrance.  The  south  jetty  is  approx¬ 
imately  1700  m  long  and  the  north  jetty  is  about  500  m 
long,  measured  from  the  shoreline.  The  ebb  shoal  is 
offset  to  the  north  and  is  of  low  relief.  The  measured 
bathymetry  does  not  include  the  water  depths  less  than 
4  m  (mean  tide  level).  Therefore  the  bathymetry  shown 
in  Fig.  18  is  a  modified  version  of  the  measured 
bathymetry  where  we  extend  the  beach  about  300  m 
to  the  shoreline  using  a  typical  beach  slope  of  1/80  in 
this  region. 
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Fig.  20.  Snapshot  of  calculated  wave-induced  currents  (77=3.5  m,  0  =  1.5,  T=  12  s). 
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A  curvilinear  grid  is  generated  using  Brackbill  and 
Saltzman’s  method  (1982)  and  is  shown  in  Fig.  19.  To 
resolve  complex  nearshore  bathymetry,  coastlines  and 
structures,  finer  resolution  is  used  near  the  coastlines, 
jetties,  and  the  inlet.  The  grid  dimension  is  140  x  121. 
The  minimum  grid  size  is  65  m  near  the  shoreline  on 
right  side  of  the  north  jetty  and  maximum  grid  size  in 
the  offshore  is  about  300  m. 

Several  test  cases  are  carried  out  using  different 
wave  conditions  chosen  according  to  the  average 
annual  storm  conditions  and  the  typical  storm  condi¬ 
tions  for  the  site.  Here  we  show  a  typical  storm  case 
with  a  wave  height  of  3.5  m  and  a  peak  wave  period  of 
12  s.  We  use  a  monochromatic  wave  in  REF/DIF-1  for 
this  case.  The  wave  direction  is  west-southwest, 
corresponding  to  7.5°  in  the  x—y  coordinates  of  the 
model. 


(a) 


Cross-shore  distance  (m)  1Q4 


Fig.  21.  Cross-shore  distributions  of  wave  height  (a),  longshore 
current  (b),  and  water  depth  (c),  along  y  =  6561  m. 


Fig.  20  shows  a  snapshot  of  the  nearshore  circula¬ 
tion  calculated  from  the  model.  Strong  longshore 
currents  develop  in  the  surfzone.  The  current  patterns 
around  the  two  jetties  are  complex  with  vortices  near 
the  tips  of  the  jetties.  The  zoom-in  figure  shows  the 
complex  patterns  of  the  current  near  the  south  jetty.  It 
should  be  noted  that  there  is  a  broad  current  flowing 
southward  near  the  tip  of  the  south  jetty  and  the  cause  of 
the  flow  needs  to  be  further  investigated.  In  the  region 
south  of  the  south  jetty,  there  are  two  peaks  of  long¬ 
shore  currents  caused  by  the  barred  beach.  Fig.  21 
shows  a  water  depth  profile  along  y  =  6 561  m  and 
corresponding  cross-shore  distributions  of  wave  height 
and  longshore  current.  The  current  magnitude  is  about 
1  m/s  which  is  a  reasonable  magnitude  of  wave-induced 
currents  under  this  wave  condition  (Gelfenbaum  et  al., 
2000). 

Though  the  numerical  results  have  not  been  com¬ 
pared  with  measurements  because  of  an  absence  of 
measurement  data  at  this  moment,  this  case  study 
illustrates  that  the  curvilinear  version  of  SHORECIRC 
has  a  potential  for  computations  in  complicated 
domains.  A  further  study  of  the  Grays  Harbor  case 
is  being  carried  out  with  considerations  of  wave- 
driven  current,  tidal  currents  and  model/data  compar¬ 
isons  and  will  be  reported  in  the  near  future. 

5.  Conclusion 

Using  a  generalized  coordinate  transformation  and 
contravariant  technique,  curvilinear  equations  of  the 
quasi-3D  nearshore  circulation  model  are  derived  in 
this  paper  based  on  the  work  of  PS99.  The  curvilin¬ 
ear  model  is  implemented  by  employing  a  high-order 
finite  difference  scheme  and  using  a  staggered  grid 
in  the  transformed  image  domain.  To  make  the 
model  applicable  to  irregular  shaped  domains  such 
as  harbors  and  tidal  inlets,  various  point  types  are 
defined  in  the  model  to  recognize  various  boundary 
conditions. 

The  numerical  convergence  is  tested,  and  very 
good  convergence  rates  with  both  space  and  time 
discretization  are  obtained.  The  model  is  then  used  in 
four  case  studies.  In  the  first  case,  a  boundary-fitted 
curvilinear  grid  is  used  in  the  simulation  of  nearshore 
circulations  generated  on  a  curved  beach  by  the 
diffracted  waves  in  the  lee  of  a  breakwater.  The 
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numerical  results  show  very  good  agreement  with 
experimental  results  by  Gourlay  (1974).  The  second 
case  shows  the  simulation  of  longshore  circulation 
around  a  conical  island.  Circulation  cells  are  found  on 
the  sides  of  the  island,  which  is  consistent  with  Mei 
and  Angelides  (1977).  The  third  case,  longshore 
current  simulation  in  an  obliquely  quadrilateral  do¬ 
main,  uses  a  non-orthogonal  Cartesian  grid  to  fit 
obliquely  oriented  wave  guides  used  in  the  LSTF 
laboratory  experiment  by  Hamilton  and  Ebersole 
(2001).  The  model/data  comparisons  of  the  longshore 
current,  wave  set-up,  as  well  as  vertical  profile  of 
cross-shore  current  show  the  accuracy  of  the  curvi¬ 
linear  model.  The  simulation  of  nearshore  circulation 
at  Grays  Harbor,  as  the  last  case,  illustrates  the 
capability  of  the  model  in  dealing  with  irregular 
geometries  in  complicated  domains.  For  this  case 
no  measured  data  is  available.  Compared  to  the 
Cartesian  version  of  the  nearshore  circulation  model 
which  is  only  used  in  rectangular  domains,  the 
present  curvilinear  version  with  the  successful  test 
cases  shows  its  versatility  to  handle  complicated 
geometries. 

Further  work  with  a  moving  shoreline  boundary 
condition  (Brocchini  et  al.,  2002)  and  an  absorbing¬ 
generating  boundary  condition  (Van  Dongeren  and 
Svendsen,  1997)  may  be  needed  to  deal  with  moving 
shoreline  problems  and  non-reflecting  seaward  bound¬ 
ary  problems.  More  practical  applications  of  the  cur¬ 
vilinear  model,  such  as  the  Grays  Harbor  case,  need  to 
be  investigated  with  field  data  verification.  As  one  of 
the  circulation  components  in  the  nearshore  commu¬ 
nity  model,  the  source  code  of  the  curvilinear  version 
of  the  Quasi-3D  nearshore  circulation  model  will  be 
available  on  line  at  http://www.coastal.udel.edu/kirby/ 
NOPP/index.html. 
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Appendix  A.  Derivations  of  equations 

A.l .  Derivations  of  Eq.  (6) 

The  variables  defined  in  Eq.  (5)  can  be  expressed 
using  the  new  basis  (a1?  a2,  k): 

u  =  i/aaa  +  wk,  V  =  Faaa  +  Wk, 

V  =  Faaa,  Vi  =  F?aa,  (39) 

uw  =  aa  +  ww  k,  u'  =  i/aaa  +  w'k.  (40) 

where  w,  and  wf  are  vertical  components  of  in¬ 
stantaneous  velocity,  short-wave  velocity  and  turbu¬ 
lence  velocity,  respectively.  V0C=V0C+  V\  and  W  is  the 
vertical  component  of  the  time  averaged  velocity. 

The  depth-integrated,  short-wave-averaged  conti¬ 
nuity  equation  may  be  written  in  a  general  form  as 

dl 

-^-  +  dwHV  =  0  (41) 

at 

where  div//  means  horizontal  divergence.  According 
to  tensor  calculus,  Eq.  (41)  can  also  be  expressed  in  a 
tensor-invariant  form: 

f-+(V«h)A  =  0  (42) 

where  ()a  is  the  covariant  derivative  of  a  first-order 
tensor.  For  an  arbitrary  vector  waaa,  for  instance,  its 
covariant  derivative  is  defined  as 
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Using  the  following  expression  of 

»  =  92x/  (gp 

*s  ~  KaKs  3x, 

Eq.  (42)  can  be  written  in  the  form  of  Eq.  (6). 
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A.2.  Derivations  of  Eq.  (9) 

The  depth-integrated,  short-wave-averaged  mo¬ 
mentum  equation  is  written  in  the  generalized  form  as 


3 

—  (Vh)  +  div// 
at 


V*Ph 


/  K 

4  —ho 


Vfdz 
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H —  div#T 
P 


-  -div//S  +  g/? gradC  -  -  ts  +  -tb  =  0  (45) 

P  P  P 


In  Eq.  (45),  S  and  T  represent  the  radiation  stress 
and  the  Reynolds’  stress,  respectively.  After  the  fol¬ 
lowing  approximation  made  in  PS99 


f  V?vfdz+  [\u«wV? +  V?Jw)dz~  P  Vfvfdz 

J-ho  j(t  J-ho 

+  vfCOQl  +  v?(OQi,  (46) 

the  simplified  form  of  the  momentum  equation  (PS99, 
Eq.  (9))  is  written  in  a  tensor-invariant  form  as 
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Usually  Sup  is  computed  from  a  wave  model  and  S 
is  then  obtained  using  Eqs.  (48) -(50).  Another  way  to 
get  the  term  Sjf  is  to  calculate  the  contravariant 
component  of  div  S  from  its  Cartesian  components 
which  is  directly  obtained  from  a  wave  model  in 
Cartesian  coordinates. 

The  term  Tff  in  Eq.  (47)  is  the  contravariant 
component  of  Reynolds’  stress  which  is  also  a  sec¬ 
ond-order  tensor.  It  can  be  calculated  using  the  rate- 
of-strain  tensor  defined  by 


D  =  I[gradF  +  (gradF)7],  (51) 

and 
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divT  =  div(AvfD),  (52) 

where  ()r  represents  the  transpose  of  a  tensor,  vt  is  the 
diffusion  coefficient.  Using  tensor  calculus,  divT  or 
Tffi  is  expressed  as 
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(53) 


In  Eq.  (47),  S"1''  is  the  contravariant  component  of 
radiation  stress  which  is  a  second-order  tensor. 
According  to  the  relationship  between  two  different 
bases,  the  relationship  between  and  Syp  can  be 
easily  obtained  as: 


where 

D^=fgy^+g^Jy)  (54) 

A.  3.  Derivation  of  Eq.  (12) 


The  general  form  of  the  horizontal  momentum 
oil  1  [c  fdxA2  oo  3x2  9X!  / 3xi  \  2\  equation  PS99  Eq.  (10)  reads 

+  (divff  (uu))ff  =  -  -  (gradp)^ 


(55) 
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where  u  =  U//+w,  ()H  means  a  horizontal  component, 
(div  (uu))//  can  be  divided  into  two  parts  below. 

(div(uu))//=div//-(u//u//)  +  (waw)z  aa  (56) 

where  a  =  1 ,2.  Extending  (naw)  z  and  using  the  tensor 
calculus 


a  \  3  uaw  o  o  „  o  3  uaw 

w)  z  =  +  r p3upw  +  r^u  uP  =  - 


3  z 


dz 


(57) 


The  derivation  above  used  the  following  values  of 
second  kind  Christoffel  symbols 

n 3=I%=0  (58) 

Then  (56)  can  be  written  as 

3  tPw 

(div(uu))^=div//(uffu//)H — - — aa  (59) 

dz 

Using  Eq.  (59),  Eq.  (55)  may  be  expressed  in  the 
tensor-invariant  form  as  shown  in  Eq.  (12) 

A.  4.  Derivation  of  Eq.  (14) 

According  to  Eq.  (59),  the  terms  involving  u'  in 
Eq.  (13)  can  be  written  in  the  general  tensor  form: 


(«,“«'/?),0  a«  +  f  («'“* ^)»«  =  (div(u,u'))//  (60) 

Reynolds  stresses  may  be  expressed  by  introducing 
the  rate-of-strain  tensor 


u'u'  =  —  v^gradF  +  (gradF)r].  (61) 

Eq.  (61)  may  be  written  in  a  tensor- in  variant  form 


as 


Li'u'  =  -vt(gyfiV*+gy*V$)  aaa^ 
3  W  3  Va 


U  g 


ya 


dL 


dz 


ka„ 


where  (a,  /?,  7)  =1,2.  The  above  derivation  used  the 
values  of  contravariant  metric: 

g3a=ga3  =  0 

g33  =  1  (63) 

and  Christoffel  symbols  (Eq.  (58)). 

Substitution  of  Eq.  (62)  into  Eq.  (60)  leads  to  the 
representation  of  ( div(u'u'))H : 


(div(u'«i'))„  =  [— v,(g^F“  +  gyxVfy)]^x 


Vt 


+  glli 


(64) 


Substituting  Eq.  (64)  into  Eq.  (13)  gives  the  tensor- 
invariant  form  (Eq.  (14)). 


Appendix  B.  Dispersive  mixing  coefficients 

Using  Eq.  (27)  and  the  expressions  of  A,  B ,  D  and 
M  in  PS99,  we  get 


Daft  — 


—  diudiph6  +  —  (di^exp  +  d\pe\f)h5 
Vf  od  do 

Y^d\oc(fip  +fip)  +  -^d\p(f\x  fif) 

^eiaei^jh4  +  — (^ia(/ij8 

e\p(f\a  +  —  (/laT, flu )  (/l  /?  +/2j?)^2 


(65) 


B. 


4  (  1  1 

—  d\pd\ah6  +  (  —diaeip-\-  —  eiocdip  )  h* 


(2  2 
+  (  —dixifip  -\~fip)  +  —dipifia  A-fi a) 

+  j^eioceip^Jh4  A-  ^-£ia(/i/l  +fip) 
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Map  —  —d\ad\ph5  +  —  (di^eip  +  dipe\a)h4 

+  ^  {d\a{flp  +flp)  +  d\p(fia  -\-f20c) 

1 

+  ^loc^2p)h  +  —  (e\a(flp  +fip) 

+  e\p(flot  +fla))h2  +  (/la  +fla)(f\p  +  flp)h 
“I-  ( Qwocd\p  Q,wfid\a)h  ~\~  (Qwoc^lp  Qw(i€\oi)h 

+  Qwotfflp  +flp)  +  Qwpiflct  +fla)  (67) 


^lj? 


]^[  —  2/z^la  j  tf?ljS  h6  +  —  JJ(/l0  -\-f2p) 

ad  )  -  L  15  ac5 


20  1 

1 

8 


^90  (  n  _ 

oa5  /  0  \a<5 

Y\  **-2hfd\u\  (f\p  +fip) 

\  a<5  / 

( n  ~h^  )eip  ^ + \  ( n  h,de\a 

\  a 5  /  -I  J  \  a<5  y 

(/ij?  +fip)h 3  +  —  ^  JJ  d\a^  h1 

.P&  \P8  J 

-f  nw«  +/2«) + 2-  ( n  )  «i 


1 

36 


A5 +77 


3  —2hjd\p 
KPS  / 


(/la  +/2a)  +  (  J^[  ~h,Se Ip  J  eio 


.P8 


l  (  II  “^1/H  (/la  +/2a)^3 


(68) 


.ps 

rl  rr2 


where  EL>  EL  and  EL  are  defined  by 
1 

"J  =  d\a^S  +  ridd\k 

ad 


(69) 


(70) 

(71) 


Appendix  C .  Staggered  grid  arrangement  and 
numerical  scheme 


A  staggered  grid  in  ^  —  £2  plane  is  employed  as 
shown  in  Fig.  22,  where  the  crosses  denote  (-points  at 
which  (is  computed,  the  circles  denote  E1 -points  at 
which  F 1  is  computed  and  the  squares  denote  F2- 
points  at  which  V2  is  computed. 

The  first-order  spatial  derivative  terms  are  discre¬ 
tized  to  fourth-order  accuracy  using  five-point  finite- 
defferencing.  For  example,  fa  may  be  discretized  at 
point  j  as 

F  =  ifj- 2  -  27 fj-x  +  27 fj  -fJ+1)/( 24A£,)  +  0(8$) 

(72) 

The  fourth-order  Adams -Bashforth- Moulton  pre¬ 
dictor-corrector  scheme  is  employed  to  perform  time 
updating.  A  sequence  of  time  instants  are  defined  by 
t=p  A t.  Level  p  refers  to  information  at  the  present, 
known  time  level.  For  the  time-derivative  equation 
written  by 


(73) 


r 

9  \ 

/  r 

9  \ 

/  r 

V 

J  / 

\  k 

J  / 

\  k 

J 

r 

y  \ 

9  \ 

k 

J  / 

\  k 

J  / 

\  k 

J 

Fig.  22.  Staggered  grid  in  -  £2  plane  (  x  £  point,  O-V1  and  □  -  V2). 
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the  predictor  step  is  the  third-order  explicit  Adams - 
Bashforth  scheme,  given  by 

fit'  =^-  +  ^[23(^y- \6{EX)  +5(E)l J2}  (74) 

After  predictor  step  we  use  the  fourth-order  Adams - 
Moulton  corrector  method: 

fir  =fij + ya  p  F)pr + i9(Er  -  mr 

+  (Etj2]  (75) 

Appendix  D.  Interpolation/extrapolation  between 
model  grids 


1 


Fig.  23.  Interpolation  triangle. 


Interpolation  or  extrapolation  is  usually  employed 
between  the  curvilinear  grid  used  for  the  circulation 
model  and  the  rectangular  grid  used  for  wave  models. 

The  linear  interpolation/extrapolation  method  is 
used  in  the  present  paper.  We  assume  two  grid  systems, 
grid-1  and  grid-2,  which  can  be  curvilinear  grids  or 
rectangular  grids.  As  shown  in  Fig.  23,  the  interpola¬ 
tion  value  at  point  A  in  grid-1  is  evaluated  by  the 
values  at  three  points,  1,  2  and  3,  of  a  triangle  in  grid-2 
which  surrounds  points.  For  extrapolation,  points  is 
out  of  the  triangle.  Four  triangle  areas  S a^7,  i.e.,  Si23, 
S12A,  S31A  and  S23a  are  calculated  using  the  following 
formula: 


xa  ya  1 
xp  yp  1 
yy  1 


(76) 


where  (xa,  ya)  represents  coordinates  of  point  1,  2,  3 
and  A.  For  interpolation,  (a,  /?,  7)  are  anti-clockwise  for 
all  the  four  triangles  and  thus  Sapy  are  positive.  For 
extrapolation,  clockwise  (a,  /?,  7)  results  in  negative 
S^py.  The  following  formula  is  used  for  both  interpo¬ 
lation  and  extrapolation: 


Fa  —  {F\S23a  +F2S$ia  +  F2S\2a)  /  S\2i  (77) 

where  F 1,  F2 ,  F3  and  FA  represent  any  converted 
variables  at  point  1,  2,  3  andv4,  respectively. 


It  should  be  mentioned  that  the  extrapolation  is 
usually  used  at  domain  boundaries  and  is  only  suitable 
for  the  case  that  domain  boundaries  are  close  to  each 
other.  For  the  cases  that  two  domain  boundaries  are 
not  close  to  each  other,  e.g.,  Gourlay’s  experiment  and 
conical  island  case,  we  make  the  interpolation  in  the 
overlapped  region  of  the  two  grids.  To  save  compu¬ 
tational  time  for  interpolation/extrapolation,  Sxpy  are 
stored  when  the  interpolation/extrapolation  subroutine 
is  first  called. 
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